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Abstract 

I present a new and reliable method to test the numerical accuracy of NLO 
calculations based on modern OPP/Generalized Unitarity techniques. A 
convenient solution to rescue most of the detected numerically inaccurate 
points is also proposed. 
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1. Introduction 

With the advent of the modern OFF [1, 2] and Generalized Unitarity [3, 
4, 5, 6, 7, 8, 9, 10] based techniques, the art of computing NLO corrections 
received a lot of attention in the last few years and several programs [11, 12, 
13, 14, 15, 16] and computations [17, 18, 19, 20, 21, 22, 23] exist, by now, 
based on this philosophy. 

While for the traditional reduction methods [24] a lot of work has been 
spent already to find ways to control the numerical accuracy of the results 
[25, 26], in the case of these new techniques the situation is still at an early 
stage. However, their potential to self detect stability problems is known 
since 2007 [11], the basic observation being that, since a reconstruction of 
a function N{q) of the would be integration momentum q is involved (the 
coefficients of which are interpreted as the coefficients of the scalar 1-loop 
functions entering the calculation), one can numerically test the accuracy of 
it by comparing N{q) and its re-constructed counterpart at a new, arbitrarily 
chosen value of q. However, the arbitrariness of the point chosen for the test 
poses serious problems, because it introduces a new, unwanted, parameter 
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upon which the check depends in an unpredictable way ^. Furthermore, not 
all the reconstructed coefficients enter into the actual computation, because, 
for example, some of them may multiply vanishing loop functions, rendering 
immaterial a possible inaccuracy in their determination. In addition, it is 
not clear how to test the rational part of the amplitude [2] . 

Nevertheless, it keeps being very tempting the idea of self detecting nu- 
merical inaccuracies, avoiding the need of additional analytic work. In this 
paper, 1 present a method to achieve this task, based on the construction of a 
reliable precision estimator working at an event by event basis. Therefore, it 
becomes possible for the user to safely set a precision threshold above which 
the inaccurate points are discarded. Furthermore, I prove that, re-fitting the 
discarded points at higher precision while keeping the computation of N{q) 
in double precision allows to re-include most of them in the original sample. 
This solution nicely factorizes the problem, in the sense that the codes (of 
the parts of the code) computing the function N{q) can be kept in double 
precision and only the fitting procedure to get the coefficients needs to be 
re-done at higher precision. 

The structure of this work is very simple: in Sections 2 and 3, I describe 
the algorithm and, in Section 4, I report on the tests I performed on the 
whole procedure. 

2. The method 

In the OPP technique the numerator N(q) of the integrand of a m-point 
amplitude is decomposed in terms of denominators Di — {q +Pi)^ — ml ^ 

m—1 m—1 m—1 m—1 

N{q) = Yl <W2;q) J] A + J] K^o^i; A 

io<ii<i2 i^io,ii,i2 io<ii i^io,ii 

m—1 m—1 



^Improvements on this technique have been recently presented in [16]. 
^In our notation, q in 4-dimensional, q n-dimensional and n-dimensional denominators 
are written as Di = (q+Pi)^ ~ ™i • 
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where, for later convenience, I have grouped all the 4-point contributions into 
a single term 

m—1 m—1 

pM(g) = J2 d{ioHi2i3;q) n (2) 

io<ii<i2<i3 i^io,ii,i2,iz 

The functions (^(io^i^2^3; c{iQiii2;q), b{iQii;q) and a{io;q) depend on the 
integration momentum q and bring information on the coefficients of the 
scalar 1-loop integrals, that are obtained by fitting N{q) at different values 
of q that nullify, in turn, the denominators. Finally, performing a global shift 
of all the masses appearing in the denominators of Eq. 1 

- f (3) 

and fitting again, allows to reconstruct also a piece of the rational terms, 
called Ri [2]. In summary, by knowing the set ^ 

(i(ioili2«3) , c(ioili2) , /^x 

b(ioii) , a(io) , Ri , 

the amplitude A is reconstructed by simply multiplying by the corresponding 
scalar 1-loop integrals [12, 27] ^. 

m—1 



.u^il<l2<«3 ° ^ ^ 

m—1 „ ^ 

+ J2 <^otit2) / d^'q 



J<ll<l2 
m—1 



m—1 „ ^ 

+ J] a{io) / + i?i . (5) 

io 



use a notation such that the cocfBcionts of the scalar 1-loop hmetions have the same 
name of the functions appearing in Eqs. 1 and 2, but without q dependence. 

^The remaining piece of the rational terms i?2 can be computed as explained in [2, 28, 
29]. 
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The key point of the method I propose in this paper is the observation 
that, if one would be able to obtain the whole set of Eq. 4 in an independent 
way, giving, as a result, a new set 

6'(zoZi), a'(io), i?;, 
an independent determination of the 1-loop amplitude would become possible 

m—1 

d'(iniTioii) I d^f, 



A' ^ J2 d'{iQHi2h) j d^'q ^ ^ ^ 

j0<«l<«2<«3 '° 
m—1 „ ^ 

io<ii<i2 *° ' 

m—1 „ 



D<»1 

m—1 



+ E «'(^o) / + K , (7) 



that could then be used to define a reliable estimator of the accuracy as 
follows ^ 

The advantage of Eq. 8, with respect to a test performed at the level of 
the function N(q), is that onZy the coefficients contributing to the amplitude 
enter into the game. Furthermore, a test on Ri becomes possible. As it will 
become clear shortly, it is convenient to differentiate the two cases where the 
coefficients of the sets in Eqs. 4 and 6 (and, a fortiori, the amphtudes in 
Eqs. 5 and 7) are computed in double or multi-precision. Then, I denote the 
double precision estimator by 

Et^^A^^ (9) 



^In an actual, numerical implementation, a small quantity e has to be included in the 
denominator of Eq. 8 to deal with the case of vanishing amplitude. 
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and its multi-precision version by 

Ei - '^7:f^' . (10) 

In the rest of this section, I illustrate how to obtain the new set of Eq. 6. 

The technique is similar to the procedure adopted to compute Under 
the shift in Eq. 3, Eq. 1 becomes 

m— 1 m—l 

N{q) = V^^\q)+ c(w2;g) H (A + 

io<i\<i2 i^io,ii,i2 
m—l m—l 

+ J2Hwq) n (Di + f) 

io<ii i^io,ii 
m—l m—l 

+ J2^('o;q)ll(Di + f), (11) 

where 

m 

with the last coefficient of f>^'^\q) independent on q 

^(2m-4)^^^_^(2m-4)^ (13) 

and where 

c{ioiii2;q) = c{ioiii2;q) + fc^^\ioiii2;q) 
b{ioii;q) = b{ioii; q) + fb^'^^ioh; q) 
a{io;q) = a{io;q). (14) 

This last equation imphes, for the l-,2 and 3-point coefficients 

c(io^i^2) = c{ioiii2) + fc'^'^\ioHi2) 
h{ioii) = b{ioii) + fb^'^\ioii) 
a{io) = a{io) . (15) 
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The constants 6^^^(io«i), c^'^\'ioHi2) and d*^^"* enter into the computation 
of Ri [2] 

m—l 
io<h<i2 

J:b^'\^o^.)(rnl+rnl-^^^^^), (16) 



327r2 

io<»i 

and can be determined with the help of Eqs. 15 and 12. 

With the knowledge of a(io), b{^oii), c(io^i^2)j ^''^■'(^o^i) a-nd c^^^(^o^l^2), 
a', b' and c' in Eq. 6 can be immediately obtained with a new mass shift 

m,' ^ m," - qf , (17) 

giving 

Ci(io^i^2) = c{ioiii2) + qlc'-^\ioiii2) 
bi{ioii) = &(«o^i) + ^i^^^H^o^i) 
ai(io) = a(io), (18) 

where I attached the subscript i to the coefficients obtained with the new 
shift. Combining Eqs. 15 and 18 gives 

a'(io) = ai{io) 

uc ■ \ Kioh) + bi{ioii) + ql ,(2)/- ■ x 
b{^o^l) = ^ ^ — h''\ioii), 

c(^o^l^2) = ^ ^ — c^>{iotit2). (19) 

As for Ri, an independent determination of b^^\ioii), c^'^\ioiii2) in Eq. 16 
also follows from the new shift 

,(2)r- ■ - A c(ioiii2) - ci(ioiii2) 
c'-'{totit2) = Z2 

^,,2, = - . (20) 

g2 _ ^2 

With the help of Eq. 18 one can now completely reconstruct the 1-, 2- and 
3-point parts of the numerator function with masses shifted according to 
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Eq. 17, namely 

m— 1 m— 1 m— 1 m— 1 

^ Ci(ioiii2;?) J]^ (Di + ql) + ^bi(ioii;q) ( A + ^i) 

m—l m— 1 

+ J]ai(zo;g)n(A + g?). (21) 

By subtracting Eq. 21 from A^(q') one determines T>'{^\q) obeying the follow- 
ing polynomial (in q\) representation (see Eq. 12) 

m 

T>r\^ = j:^}''-''di'^-'\q), (22) 

3=2 

the first coefficient of which, computed at values q = ?io,n,i2,i3 nullifying, in 
turns, all possible combinations of 4 denominators 

-^io (?io,w,j2,«3) -^hiQioMtiijia) -^12 (?io,n,i2,j3) -^isilioMjiijis) ' (^3) 

gives the desired independent determination of the box coefficients 

d'{ioiii2i3) = di°\qio,n,i^,i3) ■ (24) 
Prom the last term in Eq. 22 one obtains, instead 

^/{2m-4)^^(2m-4)^ (25) 

that, together with the coefficients in Eq. 20, gives a complete alternative 
determination of Ri 

rn—l 

R[ = L^,i'(2--4) _ ^ y c'^^\ioHi2) 

^ 967r2 327r2 . ^ . ^ " ^ 

io<ii <i2 

j:b'^^\ioH) (ml+ml- ^P--/-^' ^ , (26) 



327r2 

tO<H 

I close the Section by summarizing the procedure. One fits the numerator 
function N[q) three times; the first time with q^ = (Eq. 1) to determine 
the cut-constructible part of the amplitude, namely all of the coefficients in 
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Eq. 4; the second time with ^ 0, to compute i?i by means of Eq. 16; the 
third time with a new value of the mass shift [q^) to calculate the alternative 
set of coefficients in Eqs. 19, 20, 24 and 25, that allow to build the precision 
estimator given in Eq. 9. 

It is clear that performing the three fits by using directly the numerator 
function N{q) appearing in the l.h.s. of Eq. 1 could be computationally 
very expensive because, in practical cases, the calculation of N{q) is rather 
time consuming ®. Fortunately, after the first fit, one is allowed to use the 
reconstructed numerator function (namely the r.h.s. of Eq. 1) to determine 
both Ri and E^. The additional CPU time is then very moderate. To further 
decrease it, one can also observe that, instead of determining ci^^™"^) through 
the expansion in Eq. 12, that requires the knowledge of T>^"^\q) at (m — 2) 
different values of q'^ , one can get it by means of the following relation among 
the OPP coefficients 

m—1 m—1 m—1 m,— 1 

rf(2m-4) ^ _ J2 c^'\totit2; q) - K^o^i; q) - E ^^'H^o^i; E 

io<U<*2 io<U io<U i^io,ii 

m—1 m—1 

- ^a(zo;g)^A. (27) 
Eq. 27 is proved in Appendix A. 



3. Rescuing the inaccurate points 

Under the assumption that E^ in Eq. 9 is a good precision estimator, 
points can be rejected when E^ > -Eum, where £"11111 is a threshold value cho- 
sen by the user. In this Section, I propose, as a simple recipe to rescue the 
rejected points, to re-perform the three fits described in Section 2 at higher 
precision while keeping the computation of N{q) in double precision. The 
advantage of this recipe is that multi-precision routines need to be imple- 
mented just in the fitting program, while the code providing N{q) can be left 
untouched ^. A new test in multi-precision can then be performed on the 



assume here an unoptimized computation of N{q), performed without cashing the 
information that does not depend on q. 

^This strategy is especially relevant in the case of programs that already implement, 
internally, multi-precision routines [30], such as CutTools [11]. 
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rescued points using the multi-precision estimator of Eq. 10 and only if, even 
in this case, E;^ > E^ja the event is discarded for good (or re-computed, if 
possible, with N{q) also evaluated in multi-precision). The hope is that the 
percentage of points rejected by this second test is very limited, so that they 
can be safely ehminated from the sample. The effectiveness of this strategy 
is studied in the next Section. 

4. Testing the method 

To test the procedures described in the previous two Sections I imple- 
mented in CutTools, as a numerator function N{q) mimicking the full am- 
plitude A, one of the 120 diagrams contributing to the 77 47 scattering 
in massless QED ^. To enhance the problematic region I did not apply any 
cut on the final state particles, so that numerically unstable Phase Space 
configurations with zero Gram determinant can be freely approached. From 
a practical point of view, 1 constructed the alternative amplitude of Eq. 7 by 
keeping the 4-point coefficients and re-computing only c', b', a' and R[ with 
the help of Eqs. 19, 20 and 25. The reason is that, in practice, the derivation 
of the 4-point part of an amplitude is numerically quite stable, the bulk of 
the numerical instabilities coming from the lower-point sectors. 

Before describing in details the tests, 1 introduce, besides the estimators 
given in Eqs. 9 and Eqs. 10, a few more variables. I define two true precision 
variables as follow 

Pd=\Aa-A,\/\A,\ and = \A^ - A,\/\A,\ , (28) 

where A^ is the exact reference amplitude computed with both fits and N{q) 
in multi-precision. Pa and Pm will be used, in the following, to test the 
actual precision in the computation of A^ and Am- Furthermore, for the sake 
of comparison, I define two additional precision estimators, based on the so 
called N^N test of [11] 

E^ = \N,-N,,.,,\/\N,\ and E^ = \N, - N^^rec\/m , (29) 

where is the numerator function N{q) computed, in double precision, at a 
random value of g ^, N^^rec the same numerator, but reconstructed, in double 

*^Tliis diagram contains up to rank six 6-point functions, so it fairly represents the 
complexity of the real situations. 

picked up the point q = (1/2, -1/3, 1/4, -1/5). 



9 



1500 — 



1000 
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^r-r-i 1 1 1 1 \ 1 1 1 1 \ 1 1 1 1 h-t-^H 

-2 2 4 

Log,o(PV^f) 

Figure 1: Distribution of the ratio between the true precision = {A^ — Ae\/\Ae\ and 
two different precision estimators. The dashed histogram refers to the estimator at the 
numerator level = = |A^^ — N4^j.ec\/\Nd\, the solid one to the estimator at the 
amplitude level E^ = = \Ad - A'^\/\Ad\. 



precision, via the r.h.s. of Eq. 1, and, finally, Njn,rec is the numerator function 
reconstructed by means of a multi-precision fitting procedure. 

In Figs. 1-3, I collect the results obtained by using 3000 random, uni- 
formly distributed Phase Space points. In Fig. 1 I plot the distributions of 
the ratios Pd/E^ and Pd/E^. In the latt er case (sohd histogram), most of 
points fall within 2 orders of magnitude, indicating that E^ is expected to 
accurately estimate the true numerical precision. This is not the case for 
the estimator based on the N — N test. The long right tail in the Pd/E^ 
distribution (dashed histogram) shows that there are points for which the 
accuracy is badly overestimated by . In Fig. 2, 1 plot the distributions 
of Pd (dashed histogram) and P^ (solid histogram). It can be seen that 
keeping the computation of N{q) in double precision, while performing the 
fit in multi-precision improves the accuracy. Nevertheless, points exist for 
which Am is still not accurate enough, even if the solid plot stops order of 
magnitudes before the dashed one. In Fig. 3, 1 show the tails of the Pd 
distribution, imposing four different cuts in the value of the estimators E^ 
and E^. In the solid histograms, when E-^ > E^j^, a rescue of the point is 
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10' 




-10 -5 5 10 

Logio(Pa;) 

Figure 2: Distribution of the true precision variables = {A^ — Af.\/\Ae\ (see text). The 
dashed histograms refers to the double precision result {P^ = Pa), the solid histogram 
to the case with fitting procedure carried out in multi-precision, but numerator function 
computed in double precision [P^ = Pm)- 



performed by re-fitting the 1-loop coefficients in multi-precision, while keep- 
ing the computation of the numerator function in double precision, and, if 
also > -Eiim, the event is discarded. In the dashed histograms, the same 
procedure is apphed, but using, as estimators, and Ej^. Again, the right 
tails of the dashed histograms show that E^ and E^ are not good estima- 
tors of the numerical accuracy, while the absence of points above 10~^ in the 
case of all the solid plots, indicates that Ef and E^ are able to select the 
bad inaccurate points quite efficiently. For reader's reference I summarize, 
in tabic 1, a statistics of the number of points computed in multi-precision 
and discarded in each of the 4 cases. As a conclusion, the rescue procedure 
is able to recover most of them. 

As a final check on the goodness of the estimator at the amphtude level, 
I present, in table 2, the quantity max[Logio(Pd)] ~Logio(£'am) as a function 
of Eii^, when using the precision estimator E^. This variable measures 
the difference between the worst detected point, in an analysis like that one 
represented by the solid histograms of Fig. 3, and the chosen threshold value 
Eiim for Ef. It can be seen that, for values of Eu^ between 10~^ and 10~^, 
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^lim = 10- 



Elua = 10- 



-4-2 2 

Logio(Pd) 



-2 

Logio(Pd) 



E\im = 0.5 X 10 



-2 



^lim = 10-2 



-4 -2 



Logio(Pd) 



-4 -2 



Logio(Pd) 



Figure 3: The tails of the distributions of the true precision variable Pd, with an ad- 
ditional constraint on the value of the precision estimators ^ E^, E^ and E^ = 
\Nd — Nm,rec\/Wd\ (scc text). In the solid histograms, when E^ > E'um, a rescue of the 
point is performed by re-fitting the 1-loop coefficients in multi-precision (with numerator 
functions kept in double precision) and, if also iJ^^ > -Biim, the event is discarded. In the 
dashed histograms, the same procedure is applied, but using the estimators E^ and E;^ 
instead. 



12 





N 


Ndis 


10-4 


90 


14 


10-3 


62 


8 


.5 X 10-2 


44 


6 


10-2 


40 


6 



Table 1: The number of points computed in multi-precision (A^mp) thanks to the rescue 
procedure and the number of points discarded {N^tis) as a function of the threshold value 
Eiim- The numbers refers to the solid histograms of Fig. 3, over a total number of 3000 
events. 

overestimates the accuracy at most by 1.1 decimals and that the points 
where the overestimate is by almost 2 decimals He in the safe region of very 
small values of Enm, from which one argues that is able to detect badly 
instable Phase Space points in a reliable way. 

5. Conclusions 

I introduced a novel method to test the numerical accuracy of the NLO 
results produced by modern OPP/Generalized Unitarity techniques. The 
key ingredient is a re-computation of the 1-loop coefficients based on the 
properties of the OPP equation under a global shift of aU the masses. This 
re-computation can be performed by using the function previously recon- 
structed during the determination of the cut-constructible part of the am- 
plitude, therefore at a moderate CPU time cost. As a by-product, I also 
introduced a faster determination of one of the coefficients contributing to 
the rational piece of the amplitude. I proved, with numerical tests, the re- 
liability of the procedure and I proposed a convenient solution to rescue 
most of the detected numerically inaccurate points in a way that allows the 
computation of the integrand to remain in double precision. 
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maxfLoKmfPj)! — LoR-iaiEum) 


in-2 

-L W 


n 80 




1 1 

J. . J. 




0.74 


5 X 10-^ 


0.57 


10-^ 


1.1 


5 X 10"^ 


0.98 


10-5 


1.0 


5 X 10-6 


1.1 


10-6 


1.0 


io-« 


1.7 


10--^" 


1.8 



Table 2: The variable max[Logio(Pd)] — Logio(£^;im) as a function of Bum, when using the 
precision estimator E^. 

Appendix A. An alternative determination of d^^"*"^^ 

The last coefficient d^^"*"^) in the expansion of Eq. 12 contributes to Ri 
through Eq. 16. In this appendix, I present a novel technique to determine 
it from the other, known, coefficients of the OPP expansion. 

The starting points are Eqs. 11-14. The l.h.s. of Eq. 11 does not depend 
on g^, so that one can equate to zero the coefficients of all the powers of ^ 
appearing in the r.h.s. Prom the highest power, one obtains 

m—1 m—1 

b^^\ioii; q) + Yl «(^o; g) = , (A.l) 

io<ii io 

while the next to highest power, g(2m-4)^ gives Eq. 27. Notice that the r.h.s. 
of Eq. 27 can be computed at arbitrary values of q, allowing extra numerical 
checks. 
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